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We use the Szekeres inhomogeneous relativistic models in order to fit supernova combined data 
sets. We show that with a choice of the spatial curvature function that is guided by current observa- 
tions, the models fit the supernova data almost as well as the LCDM model without requiring a dark 
energy component. The Szekeres models were originally derived as an exact solution to Einstein's 
equations with a general metric that has no symmetries and are regarded as good candidates to 
model the true lumpy universe that we observe. The null geodesies in these models are not radial. 
The best fit model found is also consistent with the requirement of spatial flatness at CMB scales. 
The first results presented here seem to encourage further investigations of apparent acceleration 
using various inhomogeneous models and other constraints from CMB and large structure need to 
be explored next. 
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I. INTRODUCTION. 

Complementary cosmological observations have estab- 
lished that the expansion of the universe has entered a 
phase of acceleration [l], H H 0, H @, 0, . Cosmic accel- 
eration and the dark energy associated with it have been 
recognized as one of the most important and challenging 
current problems in cosmolo gy and physics, see e.g. the 
reviews E El El El El El Hi • As discussed in the lit- 
erature, cosmic acceleration can be caused by a repulsive 
dark energy, or, can be caused by some radical changes 
to gravity theory (General Relativity). However, a third 
possibility has been discussed in a number of recent pa- 
pers and that is cosmic acceleration could be an appar- 
ent effect because cosmological observables, such as the 
luminosity-distance to supernovae, are altered by inho- 
mogeneities in the universe when analyzed within exact 
inhomogeneous models. In view of the challenges pre- 
sented by the cosmic acceleration problem, all the three 
possibilities need to be thoroughly explored. 

The possibility of an apparent acceleration has been 
discussed in the literature using the inhomogeneous 
Lemaitre-Tohnan-Bondi (LTB) models, an incomplete 
list includes [H, El, HI HI H3, [gl, Hi] . In these and other- 
papers, it was shown that the LTB models can fit super- 
nova observations and the position of the first peak of the 
Cosmic Microwave Background Radiation (CMB) with- 
out any dark energy in the models. Although the LTB 
models are spherically symmetric, the results of those pa- 
pers support, as a proof of concept, that apparent accel- 
eration is a serious possibility. Thus, it is a compelling 
endeavor to explore apparent acceleration within more 
general inhomogeneous cosmological models. 

In this paper, we present a first analysis of appar- 
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ent acceleration usin g th e appealin g S zekeres inhomo- 
geneous models [H HI]; see also |27| and references 
therein for a review of the models. The models were 
originally derived by Szekeres [H [26[ as an exact solu- 
tion to Einstein's equations with a general metric that 
has no symmetries (i.e. no Killing vector fields). The 
models have been investigated analytically by several au- 
thors [H, [H, [H HI [H [H H, HE HI H| and are re- 
garded as the best exact solution candidates to represent 
the true lumpy universe we live in. For example, the 
models are put in the same classification as the observed 
lumpy universe in [43| (see the table on page 37 there). 
In reference [3l|, the authors reformulated the models 
using a coordinate system that shows that the models 
can be regarded as non-linear exact perturbations of the 
Friedmann-Lemaitre-Robertson- Walker (FLRW) models 
but where the exact non-linearity of General Relativity 
is not altered. As we discuss further, the Szekeres models 
have a flexible geometrical structure and are very likely 
to fit various cosmological data sets without the need for 
a dark energy component. 

It is perhaps worth clarifying that we are not proposing 
the Szekeres model as the true model of the universe 
but rather investigating the possibility of an apparent 
acceleration in an inhomogeneous universe. 

In this scenario, apparent acceleration is due to the 
fact that we happen to live in one of the many under- 
dense regions of the universe. The physical reason for 
such an observation is that inside the under-dense re- 
gion, there is less matter to slow down the expansion 
and therefore the expansion rate in the under-dense re- 
gion is larger than what it is in an over-dense region 
or compared to the overall averaged rate of expansion. 
Some limits have been put on the dispersion in this ex- 
pansion, see for example [39J, but this was contested, see 
for example [4fj, |4l|, |42| . So the difference between the 
dynamics inside the under-dense region compared to out- 
side of it could give the apparent effect of an accelerating 
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expansion. Furthermore, in these models we do not have 
to impose that we are located close to the center of the 
under-dense region because the definition of a center it- 
self depends on a spherically symmetric geometry which 
our models do not suffer from. 

Interestingly, the possibility of apparent acceleration 
connects to the avera gin g problem in relativity and cos- 
mology, e.g. [H, l45l l4(jl]. The problem states that the 
operation of smoothing inhomogeneities in the universe 
from small scales of distance to large scales is an oper- 
ation that does not commute with the operation of ap- 
plying the Einstein equations. In other words, the or- 
der in which the two operations are applied is important 
and leads to different results due to the non-linearity of 
Einstein's equations. Because of this non-commutation 
problem, it becomes clear why using observables in inho- 
mogeneous models, such as the Szckcrcs models, is dif- 
ferent from using FLRW plus perturbations and why it 
can lead to different interpretations. The problem was 
pointed out well before the cosmic acceleration problem 
and has now re-emerged because of dark energy and pos- 
sible links to inhomogeneous models; for a recent review 
see (47j and references therein. In our work here, non- 
linearities affect observables and lead to apparent accel- 
eration, and this is different from other ideas where au- 
thors tried to use non-linearities in order to generate a 
true acceleration; see also [47[ for a review. 

Finally, it is worth mentioning that the possibil- 
ity of apparent acceleration also addresses the accel- 
eration/dark energy coincidence problem, e.g. reviews 
[I El El El ES EI El ■ This problem is stated in liter- 
ature as the "why now?" problem, i.e. why cosmic accel- 
eration is taking place during the current epoch of cosmic 
evolution and only after initial structures had formed? 
As pointed out in previous literature, e.g. [47| and refer- 
ences therein, in the apparent acceleration scenario, the 
effect is taking place at the current epoch because of the 
onset of non-linear cosmic structures in the universe. 

The usual difficulty in comparing inhomogeneous mod- 
els with observations is that it is not possible, in general, 
to analytically solve the null geodesic equation in these 
models and derive observable functions. This equation is 
easily solved in the FLRW model, but that is not the case 
in more complex models. In the Szekeres models, unlike 
the LTB and the FLRW models, the null geodesies are 
not radial so the full set of equations need to be inte- 
grated in order to consider observations. In the present 
work, we explore a numerical integration in order to over- 
come the problem. 

In this first of a series of papers, we launch a program 
where we plan to use relativistic cosmological models that 
are more complex than the FLRW models (focusing on 
the Szekeres and Szekeres-Szafron models) and compare 
them to currently available and future cosmological data. 
We show in this paper that the Szekeres models success- 
fully fit supernova data and provide preliminary support 
to the possibility of apparent acceleration. A number of 
points need to be further investigated and are subject to 



ongoing and future work by the authors. 

II. THE SZEKERES COSMOLOGICAL MODELS 

The models are described in some detail in [2?], H3| 
and we give here a brief overview only. We use the set of 
coordinates used by [H, HU and by [34[ . The spacetime 
metric reads 



ds 2 



-dt 



2 , (# dr 2 + R 2d P 2 + d? 



e - k(r) 



E 2 



(1) 



where R = R(t, r) is the area distance and ' is used for 
d/dr. The function k(r) is related to the energy per unit 
mass and determines the curvature of the spatial sections 
t =constant. This function divides the models into sub- 
cases: hyperbolic (fc(r) < 0), parabolic (fc(r) = 0), and 
elliptic (k(r) > 0). The function E = E(r,p,q) is given 
by 

wt , (P 2 + g 2 ) P(r) Q(r) , f0 . 

and the functions P(r), S(r),Q(r) and C(r) satisfy the 
relation 



P\r) Q 2 (r) S(r) 
C[r) - 2S(r) + 2S(r) + 2 £ ' 



(3) 



but are otherwise arbitrary. The geometrical function 
E(r,p,q) and the geometrical constant e = 0, +l,or —I 
define further sub-cases of the Szekeres models and con- 
trol the mapping of their various hyper-surfaces [34j . 
The Einstein field equations with a dust source are 



*<*■'>- Mr! -*M 

R{t,r) 



and 



8wp(t,r,p,q) = 



2(M'(r)-3A/(r)f) 
R 2 (R'-R§) 



(4) 



(5) 



where the function M(r) represents the total active grav- 
itational mass in the case e = +1 [H,|3^{. The evolution 
of R(t, r) depends on k(r) and is given by: 
Hyperbolic case: k(r) < 



R(t,r) = M }P (cosh 77 - I) 



t-t B {r) 



-k(r) 



M(r) 



-fc(r) 3 / 2 



(sinh?7 — 77) 



Parabolic case: k{r) = 



R(t,r) = M{r) 



(6) 



(7) 



(8) 
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t-t B (r) =M{r) 



Elliptic case: k(r) > 



/ x Mir), 
R(t,r) - 1- cos 77 

k[r) 



(9) 



(10) 
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_2( R \( E < rE >q 



~ E 2- 



EE, 



R r — RE r /E 



-(-) 

K E 2 > \R. 



1 - k 



- RE }r /E 



E r-E p — EE 



pr 



R r — RE r I E 



WW 



( dX )[ dX' 



dX' 



= 



t-t B (r) 



M(r) 
(fc(r)) 3 / 2 



(77 — sin r;) 



(11) 



where t B (r) is an arbitrary function of r and represents 
the Big Bang time. For a simultaneous Big Bang, we 
choose t B — constant = just as in an FLRW model. 
It can be seen that there are several sub-cases of the 
Szekeres models and that a given model is specified by 
six functions that can be reduced to five by using the 
coordinate freedom in r. The model we chose for our 
analysis is specified in section IV. 



III. THE NULL GEODESIC EQUATIONS IN 
THE SZEKERES MODELS 



The null geodesic equation describes the motion of 
light rays arriving to us from astronomical objects, and 
it is necessary to solve it in order to derive observables, 
such as the luminosity-distance to supernovae or the an- 
gular diameter distance to the CMB last scattering sur- 
face. This equation is easily solved in the FLRW but not 
in the Szekeres models, and here we employ a numerical 
approach to the problem. 

Furthermore, it is worth noting that unlike in the LTB 
geometry, there are no radial null geodesies in the Szek- 
eres models as, for example, was discussed in [48l l49ll50|. 
This is due to the dependence of these geodesies on the 
p and q coordinates and one must integrate the full set 
of geodesic equations as we proceed below. 

First, we apply the null geodesic equation 



k a - b k b = 



(12) 



using the Szekeres metric (JXJ) and where k a is a null vec- 
tor (k a k a = 0) tangent to the geodesies. The resulting 
equations for the coordinates {£, r, p, q} are written be- 
low. These were previously given in for example (48l . |49| 
and we rederived them here using the computer algebra 
system GRTensorll [Hj]: 



t R t r — R fE ? 
dX? + 
RRt 



E 2 



1 - k 

dp, 2 , / da 



R. r — R- 



E 



K dX' 
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fr_ 

dX 2 



R, 
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dt dr 

^dx'^dx' 

r /E + R{E, r /E) 2 



d2 P 1 nl R ^M dt M d P- 

dX 2+ [ R ,[ dX ,[ dX' 



(15) 



1 R r — RE r /E 
R 



1-fc 
E 



(E^ r E_p — EE j 



+ V R E ) [ dX >[ dX> [ E )( dX> 



(16) 



1 R r — RE r /E 
R~ 

'R, 



1 - k 

E 



(^E r Eq — EE 



qr , 



! 2 



E. 
~E 

d q,2 



>dX' 





where A is a parameter. 
k a k a = 0) reads 



Next, the first integral (i.e. 



( dx' 



+ 



(R r - RE, r /E) 2 dr 



R 



, dq , 



k 



K dX' E^^dX 1 K dX> ' 

(17) 

In order to study the propagation of null rays in the 
Szekeres models, one need to integrate numerically the 
set of equations (fT3"| - (fll))) subject to equation (fT7|) . 

Before to do that, let's relate these equations to the 
redshift. For that, we consider two consecutive wave- 
crests of a light ray traveling to the observer. One is 
described by the first integral at t(X) and the second is 
at t(X) + T(A), where T is, therefore, the period of the 
light. Thus we have: 



[ d\' 



(R,T 



-RE, r /Ef dr a R*_ (( dp,2 ,dq^ 
K dX> + E* [( dX> +[ dX> ' 



d{t + T) s 
dX ' 



(R(t+T,r). r -R(t+T,r)E, r /E) 2 



1-k 



I //dp 



(18) 
(f ) 2 (19) 



) 2 + (^) 2 ) 



We can now use the first integral to see how the period 
changes as the light propagates. We will assume that the 
period of the light, T, is much smaller than the time over 

(14) 

which the light propagates, see for example [50]. This 
allows us to Taylor expand to first order 



R r — RE r J E 



R(t + T,r) = R{t, r) + R{t, r)T 



(20) 
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and 



R'(t + T,r) = R'(t, r) + R'(t, r)T 



(21) 



substituting into (| 18[) and (1201) , and subtracting (|T8|) from 
([20]) we get: 



d\ ~dX 



R' R+RR( j|-) 2 — (_/?,' R+RR' 
T^fc 



f?((^) 2 + (l!) 2 )) 



)(^r) 2 (22) 



Now, redshift is related to T by 



l + z(X e ) 



T(Xo) 

T(Xe) 



(23) 



where the subscript e is at emission and o is at observa- 
tion. We take the derivative with respect to A e to obtain: 



dz 

dX e 

We can now substitute 
scripts to get: 



dT{X e ) 1 + z(X e ) 



(24) 



dX e T(A e ) 

into (f2"3"]) and drop the sub- 



dz 
dX 



-(1 



Z)( ^ )(( «2W 



-) 2 -(R,'R,+RR')i 



l-k 

-#((^) 2 + (H) 2 )) 



: )(ix) 2 



(25) 



or more simply: 
d(ln(l + z)) 



dX 



1 R' R+RR( ) 2 - (R' R+RR' ) jfc w dr v 2 
ifu 1-fc AdA- 1 



l?((^) 2 + (^) 2 )) 



(26) 



Now, equations (|T3ll - (|16ll constitute a system of 
second-order ordinary differential equations (ODEs) for 
the functions {t(X),r(X),p(X),q(X)} where the coeffi- 
cients are composed of the metric functions evaluated 
on the null cone using the Field equations and the model 
specifications. Equation (|26|) relates the system to the 
redshift. 

We integrate the system using a fourth-order Runge- 
Kutta algorithm with adaptive step size 53]. Our nu- 
merical code iterates between calls to evaluate the Field 
equations on the null cone and calls to integrate the 
ODEs. We also coded a 3-dimensional interpolator [53| 
in order to select geodesies with starting points matching 
the coordinates of the supernovae of interest as given in 
[55l ]. We implemented the Runge-Kutta code [Hj| with 
the function vectors 



and 



dy 
dl 



{ 



, dt dr dp dq. 



dt dr dp dq d 2 t d 2 r d 2 p d 2 q 

li , ll , ll , li^ , li 2, ll 2, ~di 2 



} 



(27) 



(28) 



so that the system of 4 second-order ODEs is reduced 
into a system of 8 first-order ODEs. We also use the 



parameter I = ln(l + z) instead of A. Further detail 
about the numerical integration and the 3D-interpolator 
will be presented in a separate follow up paper [52| ■ 

The next step of our numerical integration is to eval- 
uate numerically the area distance and the luminosity 
distance. For a general definition of the observer area 
distance, ro, we go back to earlier work by [45l.l56| where 
for a general metric in spherical coordinates, ro is given 
by 



;sin 2 6i = det^J 



(29) 



where H l; is the 9-<f> block of the metric. For that, it 
is useful to rewrite the Szekeres metric in spherical co- 
ordinates. For this purpose we consider the following 
coordinate transformation: 



V = cot ( - ) cos(0), 



c.»! [ - ] sm(4>), 



(30) 



that yields the new line element 



ds 2 



-dt 2 



e — k{r) 



-dr 2 



\dd 2 +sin 2 (6>)d</> 2 ) 
(31) 

where the function E = E(l — cos(9)) does not alter other 
parts of the metric. 

Now, we note from the metric (|3ip that there are no 
off-diagonal components of the {8, (j)} part of the metric, 
i-e ggc/) = 0, and since we specialize here our analysis to 
the case of Szekeres models with e = +1, then, by P9"j) , 
the area distance can be described by [5J| 



ro = 



R 
E' 



(32) 



which indeed depends on the three coordinates r, 9 and (f> 
on the null cone. It follows that the luminosity distance 
is also a function of the three coordinates and is given by 



d L (z) = (l + z) 2 §- 



(33) 



We integrate numerically d^ along with the system of 
ODEs (fT3| - (fl6|) above. Then, as usual, we evaluate the 
distance modulus to supernova as 



m(z)-M = 51og(d L ) + 25 



(34) 



and we use the supernova magnitudes and uncertainties 
from the combined data sets of [57|, [53, [59(j and we trans- 
form the coordinates from the supernova database (55|. 
The findings are presented in the next section. 



IV. RESULTS AND DISCUSSION 

We discuss here our results from the numerical inte- 
gration and comparison to combined supernova data sets 
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From SN data set used by Davis et al. 2007 
Szekeres model 
LCDM model 
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FIG. 1: Supernova fits for the Szekeres model (green crosses) 
and LCDM (blue curve) models. The data is 94 Supernova 
(up to 1 + z = 1.449) from Davis et al 2007, Wood-Vasey et 
al 2007, and Riess et al 2007 [13, H, H| . The Szekeres model 
fits the data with a \ 2 — 112. This is close to the \ 2 — 105 
of the LCDM concordance model. Because of the possible 
systematic uncertainties in the supernova data, it is not clear 
that the difference between the two \ 2 i s significant. Fur- 
thermore, other Szekeres models remain to be explored. The 
Szekeres model used is also consistent with the requirement 
of spatial flatness at CMB scales. 



of [57|, l58|, l59(. We specify our Szekeres model by us- 
ing some guidance from current observations of matter 
abundances, CMB observations, and also some previous 
works on inhomogeneous models [21L l36l |38| . 

• We consider the case of Szekeres models with e = 
+1. This case includes the three sub-cases of hy- 
perbolic, elliptic, and parabolic models. 

• We set i_e(r) =0 for a simultaneous Big Bang like 
in the FLRW models. 

• In the usual way, e.g. [34j . we use the coordinate 
freedom in r and set M(r) — (sinh(r)) 3 . 

• In a similar way to the work done on the LTB mod- 
els by [IJl, we set here, for the Szekeres models 
k{r) = k(r)r 2 with 



jfe(r) = 



-1 



1 + C 2 r 2 



(35) 



So at large r, the spatial curvature goes to zero as 
indicated by CMB observations [{J [1| and at very 
small r curvature goes negative in accord with local 
observations of matter abundances with a density 
about one quarter of the critical density. The con- 
stant C will be set by best fit model to the data. 

For {S,P,Q}, we use model-1 of H with 
{140, 10, —113 ln(l+r)} that was considered within 
the context of structure formation using Szekeres 



and also within a general discussion of cosmological 
applications of the models [37| • Other sub- models 
will be explored in forthcoming investigations. 

Our fitings are presented in FIG. 1. We find that 
the Szekeres model that we considered fits the supernova 
data competitively up to a redshift 1 + z = 1.449. The 
best fit model has a value C — 0.80 with a x 2 — H2, 
and this is found to be close to the x 2 — 105 of the fiat 
LCDM model with Sl m = 0.27 and Q\ = 0.73 (see also 
[57[ | for the LCDM model). In view of the possible sys- 
tematic uncertainties involved in the supernova data, it is 
not clear that the difference between the two x 2 1S signif- 
icant. Also, the model that we used here was derived and 
discussed from fits that modeled large scale structure (36[ 
(clusters and voids), so while the model is found to pro- 
vide a good fit to the supernova data with (1 + z < 1.45), 
one needs to take into account the transition to larger 
scales. Indeed, it is well known that the Szekeres models 
average to almost Friedmann models at very large scale of 
distances [IE [H, [H, H3, HO, HH, H3, Hl| but further work 
is needed in our numerical approach in order to model 
such a transition and include it into the fit. Neverthe- 
less, as shown in equation ([33]) . the constant C controls 
spatial curvature and our model is consistent with the 
requirement of spatial flatness at the CMB scales. 

Now, while the LCDM model requires a cosmological 
constant in order to fit the data, the Szekeres model 
is used without any dark energy component. As indi- 
cated by the function k(r) in (|35[) . the interpretation of 
this result is that we may happen to live in one of the 
many under-dense regions of the universe as described by 
a Szekeres inhomogeneous cosmological model and this 
would leads to an apparent acceleration. 

However, a number of points need to be investigated in 
order to fully explore this possibility and work is ongoing 
by the authors. These include: 

• Exploration of other sub-models and ansatz for the 
functions S, P, Q and the curvature function k(r) of the 
model including the transition from Szekeres models to 
an almost FLRW models 

• Comparison to CMB observables. We expect that the 
position of the first peak of the CMB spectrum can be 
easily met by the Szekeres model that we considered. 
Most importantly, the models show promise to fit the 
full CMB spectra because of their flexible geometrical 
structure. Another constraint that need to be studied is 
the amplitude of primordial fluctuations in the gravita- 
tional field. Here, one needs to use the Szekeres-Szafron 
models [H, [HJ that include pressure as well. 

• Comparison of the models to observables of large- 
scale structure growth, such as gravitational lensing and 
galaxy clustering 

In conclusion, this first work of a series shows that the 
Szekeres inhomogeneous models provide a good fit to the 
supernova data and seems to support apparent acceler- 
ation but further investigations using these models are 
required in order to thoroughly explore this possibility. 
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